/*
Copyright (c) 2012 Eric Lasota
Copyright (c) 1997 Jung woo-jae

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Library General Public
    License as published by the Free Software Foundation; either
    version 2 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Library General Public License for more details.

    You should have received a copy of the GNU Library General Public
    License along with this library; if not, write to the Free
    Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/


#define _USE_MATH_DEFINES

#include <math.h>
#include <string.h>
#include <stdio.h>
#include "config.h"
#ifdef MP2DEC_ENABLE_SSE2
#include <emmintrin.h>
#endif

#include "mp2dec_internal.h"

#define PVERBOSE 0

#define MP2DEC_ALIGN_OFFSET(ptr) (MP2DEC_SIMD_SIZE - ((ptr) - static_cast<const char *>(NULL)) % MP2DEC_SIMD_SIZE)

#define MP2DEC_BIT_MASK(offset, length)  ((~((-1) << (length))) << (offset))
#define MP2DEC_STREAM_BYTE(n) (stream_reverse[n])
#define MP2DEC_STREAM_BIT(byte, bitOffs) ((byte & (1 << (7 - (bitOffs)))) >> (7 - (bitOffs)))

#define MAXSUBBAND	32
#define MAXCHANNEL  2
#define SCALEBLOCK  12
#define CALCBUFFERSIZE	512
#define MAXHISTORY		16
#define SCALE          32768
#define MAXSCALE       (SCALE-1)
#define MINSCALE       (-SCALE)

static mp2dec_uint8 stream_reverse[256];

static MP2DEC_SIMD_ALIGN FixedReal filter[16][32]=
{
	{
		 0.000000000,-0.000015259,-0.000015259,-0.000015259,-0.000015259,-0.000015259,-0.000015259,-0.000030518,
		-0.000030518,-0.000030518,-0.000030518,-0.000045776,-0.000045776,-0.000061035,-0.000061035,-0.000076294,
		-0.000076294,-0.000091553,-0.000106812,-0.000106812,-0.000122070,-0.000137329,-0.000152588,-0.000167847,
		-0.000198364,-0.000213623,-0.000244141,-0.000259399,-0.000289917,-0.000320435,-0.000366211,-0.000396729,
	},
	
	{
		-0.000442505,-0.000473022,-0.000534058,-0.000579834,-0.000625610,-0.000686646,-0.000747681,-0.000808716,
		-0.000885010,-0.000961304,-0.001037598,-0.001113892,-0.001205444,-0.001296997,-0.001388550,-0.001480103,
		-0.001586914,-0.001693726,-0.001785278,-0.001907349,-0.002014160,-0.002120972,-0.002243042,-0.002349854,
		-0.002456665,-0.002578735,-0.002685547,-0.002792358,-0.002899170,-0.002990723,-0.003082275,-0.003173828,
	},
	
	{
		 0.003250122, 0.003326416, 0.003387451, 0.003433228, 0.003463745, 0.003479004, 0.003479004, 0.003463745,
		 0.003417969, 0.003372192, 0.003280640, 0.003173828, 0.003051758, 0.002883911, 0.002700806, 0.002487183,
		 0.002227783, 0.001937866, 0.001617432, 0.001266479, 0.000869751, 0.000442505,-0.000030518,-0.000549316,
		-0.001098633,-0.001693726,-0.002334595,-0.003005981,-0.003723145,-0.004486084,-0.005294800,-0.006118774,
	},
	
	{
		-0.007003784,-0.007919312,-0.008865356,-0.009841919,-0.010848999,-0.011886597,-0.012939453,-0.014022827,
		-0.015121460,-0.016235352,-0.017349243,-0.018463135,-0.019577026,-0.020690918,-0.021789551,-0.022857666,
		-0.023910522,-0.024932861,-0.025909424,-0.026840210,-0.027725220,-0.028533936,-0.029281616,-0.029937744,
		-0.030532837,-0.031005859,-0.031387329,-0.031661987,-0.031814575,-0.031845093,-0.031738281,-0.031478882,
	},
	
	{
		 0.031082153, 0.030517578, 0.029785156, 0.028884888, 0.027801514, 0.026535034, 0.025085449, 0.023422241,
		 0.021575928, 0.019531250, 0.017257690, 0.014801025, 0.012115479, 0.009231567, 0.006134033, 0.002822876,
		-0.000686646,-0.004394531,-0.008316040,-0.012420654,-0.016708374,-0.021179199,-0.025817871,-0.030609131,
		-0.035552979,-0.040634155,-0.045837402,-0.051132202,-0.056533813,-0.061996460,-0.067520142,-0.073059082,
	},
	
	{
		-0.078628540,-0.084182739,-0.089706421,-0.095169067,-0.100540161,-0.105819702,-0.110946655,-0.115921021,
		-0.120697021,-0.125259399,-0.129562378,-0.133590698,-0.137298584,-0.140670776,-0.143676758,-0.146255493,
		-0.148422241,-0.150115967,-0.151306152,-0.151962280,-0.152069092,-0.151596069,-0.150497437,-0.148773193,
		-0.146362305,-0.143264771,-0.139450073,-0.134887695,-0.129577637,-0.123474121,-0.116577148,-0.108856201,
	},
	
	{
		 0.100311279, 0.090927124, 0.080688477, 0.069595337, 0.057617188, 0.044784546, 0.031082153, 0.016510010,
		 0.001068115,-0.015228271,-0.032379150,-0.050354004,-0.069168091,-0.088775635,-0.109161377,-0.130310059,
		-0.152206421,-0.174789429,-0.198059082,-0.221984863,-0.246505737,-0.271591187,-0.297210693,-0.323318481,
		-0.349868774,-0.376800537,-0.404083252,-0.431655884,-0.459472656,-0.487472534,-0.515609741,-0.543823242,
	},
	
	{
		-0.572036743,-0.600219727,-0.628295898,-0.656219482,-0.683914185,-0.711318970,-0.738372803,-0.765029907,
		-0.791213989,-0.816864014,-0.841949463,-0.866363525,-0.890090942,-0.913055420,-0.935195923,-0.956481934,
		-0.976852417,-0.996246338,-1.014617920,-1.031936646,-1.048156738,-1.063217163,-1.077117920,-1.089782715,
		-1.101211548,-1.111373901,-1.120223999,-1.127746582,-1.133926392,-1.138763428,-1.142211914,-1.144287109,
	},
	
	{
		 1.144989014, 1.144287109, 1.142211914, 1.138763428, 1.133926392, 1.127746582, 1.120223999, 1.111373901,
		 1.101211548, 1.089782715, 1.077117920, 1.063217163, 1.048156738, 1.031936646, 1.014617920, 0.996246338,
		 0.976852417, 0.956481934, 0.935195923, 0.913055420, 0.890090942, 0.866363525, 0.841949463, 0.816864014,
		 0.791213989, 0.765029907, 0.738372803, 0.711318970, 0.683914185, 0.656219482, 0.628295898, 0.600219727,
	},
	
	{
		 0.572036743, 0.543823242, 0.515609741, 0.487472534, 0.459472656, 0.431655884, 0.404083252, 0.376800537,
		 0.349868774, 0.323318481, 0.297210693, 0.271591187, 0.246505737, 0.221984863, 0.198059082, 0.174789429,
		 0.152206421, 0.130310059, 0.109161377, 0.088775635, 0.069168091, 0.050354004, 0.032379150, 0.015228271,
		-0.001068115,-0.016510010,-0.031082153,-0.044784546,-0.057617188,-0.069595337,-0.080688477,-0.090927124,
	},
	
	{
		 0.100311279, 0.108856201, 0.116577148, 0.123474121, 0.129577637, 0.134887695, 0.139450073, 0.143264771,
		 0.146362305, 0.148773193, 0.150497437, 0.151596069, 0.152069092, 0.151962280, 0.151306152, 0.150115967,
		 0.148422241, 0.146255493, 0.143676758, 0.140670776, 0.137298584, 0.133590698, 0.129562378, 0.125259399,
		 0.120697021, 0.115921021, 0.110946655, 0.105819702, 0.100540161, 0.095169067, 0.089706421, 0.084182739,
	},
	
	{
		 0.078628540, 0.073059082, 0.067520142, 0.061996460, 0.056533813, 0.051132202, 0.045837402, 0.040634155,
		 0.035552979, 0.030609131, 0.025817871, 0.021179199, 0.016708374, 0.012420654, 0.008316040, 0.004394531,
		 0.000686646,-0.002822876,-0.006134033,-0.009231567,-0.012115479,-0.014801025,-0.017257690,-0.019531250,
		-0.021575928,-0.023422241,-0.025085449,-0.026535034,-0.027801514,-0.028884888,-0.029785156,-0.030517578,
	},
	
	{
		 0.031082153, 0.031478882, 0.031738281, 0.031845093, 0.031814575, 0.031661987, 0.031387329, 0.031005859,
		 0.030532837, 0.029937744, 0.029281616, 0.028533936, 0.027725220, 0.026840210, 0.025909424, 0.024932861,
		 0.023910522, 0.022857666, 0.021789551, 0.020690918, 0.019577026, 0.018463135, 0.017349243, 0.016235352,
		 0.015121460, 0.014022827, 0.012939453, 0.011886597, 0.010848999, 0.009841919, 0.008865356, 0.007919312,
	},
	
	{
		 0.007003784, 0.006118774, 0.005294800, 0.004486084, 0.003723145, 0.003005981, 0.002334595, 0.001693726, 
		 0.001098633, 0.000549316, 0.000030518,-0.000442505,-0.000869751,-0.001266479,-0.001617432,-0.001937866,
		-0.002227783,-0.002487183,-0.002700806,-0.002883911,-0.003051758,-0.003173828,-0.003280640,-0.003372192,
		-0.003417969,-0.003463745,-0.003479004,-0.003479004,-0.003463745,-0.003433228,-0.003387451,-0.003326416,
	},
	
	{
		 0.003250122, 0.003173828, 0.003082275, 0.002990723, 0.002899170, 0.002792358, 0.002685547, 0.002578735,
		 0.002456665, 0.002349854, 0.002243042, 0.002120972, 0.002014160, 0.001907349, 0.001785278, 0.001693726,
		 0.001586914, 0.001480103, 0.001388550, 0.001296997, 0.001205444, 0.001113892, 0.001037598, 0.000961304,
		 0.000885010, 0.000808716, 0.000747681, 0.000686646, 0.000625610, 0.000579834, 0.000534058, 0.000473022,
	},
	
	{
		 0.000442505, 0.000396729, 0.000366211, 0.000320435, 0.000289917, 0.000259399, 0.000244141, 0.000213623,
		 0.000198364, 0.000167847, 0.000152588, 0.000137329, 0.000122070, 0.000106812, 0.000106812, 0.000091553,
		 0.000076294, 0.000076294, 0.000061035, 0.000061035, 0.000045776, 0.000045776, 0.000030518, 0.000030518,
		 0.000030518, 0.000030518, 0.000015259, 0.000015259, 0.000015259, 0.000015259, 0.000015259, 0.000015259,
	}
};

static FixedReal scalefactorstable[64];

static FixedReal group5bits[27*3];
static FixedReal group7bits[125*3];
static FixedReal group10bits[729*3];


static const FixedReal *grouptableA[16] =
{ 0,group5bits,group7bits,group10bits,0,0,0,0,0,0,0,0,0,0,0,0};
static const FixedReal *grouptableB1[16] =
{ 0,group5bits,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
static const FixedReal *grouptableB234[16] =
{ 0,group5bits,group7bits,0,group10bits,0,0,0,0,0,0,0,0,0,0,0};

static const int codelengthtableA[16] =
{ 0, 5, 7, 10, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 };
static const int codelengthtableB1[16] =
{ 0, 5, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 };
static const int codelengthtableB2[16] =
{ 0, 5, 7, 3, 10, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 16 };
static const int codelengthtableB3[8] = { 0, 5, 7, 3, 10, 4, 5, 16 };
static const int codelengthtableB4[4] = { 0, 5, 7, 16 };


static const FixedReal factortableA[16] =
{ 0.0,        1.0/2.0,    1.0/4.0,    1.0/8.0,
  1.0/8.0,    1.0/16.0,   1.0/32.0,   1.0/64.0,
  1.0/128.0,  1.0/256.0,  1.0/512.0,  1.0/1024.0,
  1.0/2048.0, 1.0/4096.0, 1.0/8192.0, 1.0/16384.0 };
static const FixedReal factortableB1[16] =
{ 0.0,        1.0/2.0,    1.0/4.0,     1.0/8.0,
  1.0/16.0,   1.0/32.0,   1.0/64.0,    1.0/128.0,
  1.0/256.0,  1.0/512.0,  1.0/1024.0,  1.0/2048.0,
  1.0/4096.0, 1.0/8192.0, 1.0/16384.0, 1.0/32768.0 };
static const FixedReal factortableB2[16] =
{ 0.0,        1.0/2.0,    1.0/4.0,    1.0/4.0,
  1.0/8.0,    1.0/8.0,    1.0/16.0,   1.0/32.0,
  1.0/64.0,   1.0/128.0,  1.0/256.0,  1.0/512.0,
  1.0/1024.0, 1.0/2048.0, 1.0/4096.0, 1.0/32768.0 };
static const FixedReal factortableB3[8] =
{ 0.0, 1.0/2.0, 1.0/4.0, 1.0/4.0, 1.0/8.0, 1.0/8.0, 1.0/16.0, 1.0/32768.0 };
static const FixedReal factortableB4[4] = { 0.0, 1.0/2.0, 1.0/4.0, 1.0/32768.0 };


static const FixedReal ctableA[16]=
{           0.0, 1.33333333333, 1.60000000000, 1.77777777777,
  1.06666666666, 1.03225806452, 1.01587301587, 1.00787401575,
  1.00392156863, 1.00195694716, 1.00097751711, 1.00048851979,
  1.00024420024, 1.00012208522, 1.00006103888, 1.00003051851};
static const FixedReal ctableB1[16]=
{           0.0, 1.33333333333, 1.14285714286, 1.06666666666,
  1.03225806452, 1.01587301587, 1.00787401575, 1.00392156863,
  1.00195694716, 1.00097751711, 1.00048851979, 1.00024420024,
  1.00012208522, 1.00006103888, 1.00003051851, 1.00001525902};
static const FixedReal ctableB2[16] =
{           0.0, 1.33333333333, 1.60000000000, 1.14285714286,
  1.77777777777, 1.06666666666, 1.03225806452, 1.01587301587,
  1.00787401575, 1.00392156863, 1.00195694716, 1.00097751711,
  1.00048851979, 1.00024420024, 1.00012208522, 1.00001525902};
static const FixedReal ctableB3[8] =
{           0.0, 1.33333333333, 1.60000000000, 1.14285714286,
  1.77777777777, 1.06666666666, 1.03225806452, 1.00001525902 };
static const FixedReal ctableB4[4] =
{ 0.0, 1.33333333333, 1.60000000000, 1.00001525902 };


static const FixedReal dtableA[16]=
{           0.0, 0.50000000000, 0.50000000000, 0.50000000000,
  0.12500000000, 0.06250000000, 0.03125000000, 0.01562500000,
  0.00781250000, 0.00390625000, 0.00195312500, 0.00097656250,
  0.00048828125, 0.00024414063, 0.00012207031, 0.00006103516};

static const FixedReal dtableB1[16]=
{           0.0, 0.50000000000, 0.25000000000, 0.12500000000,
  0.06250000000, 0.03125000000, 0.01562500000, 0.00781250000,
  0.00390625000, 0.00195312500, 0.00097656250, 0.00048828125,
  0.00024414063, 0.00012207031, 0.00006103516, 0.00003051758};

static const FixedReal dtableB2[16]=
{ 0.0,           0.50000000000, 0.50000000000, 0.25000000000,
  0.50000000000, 0.12500000000, 0.06250000000, 0.03125000000,
  0.01562500000, 0.00781250000, 0.00390625000, 0.00195312500,
  0.00097656250, 0.00048828125, 0.00024414063, 0.00003051758};

static const FixedReal dtableB3[8]=
{           0.0, 0.50000000000, 0.50000000000, 0.25000000000,
  0.50000000000, 0.12500000000, 0.06250000000, 0.00003051758};

static const FixedReal dtableB4[4]=
{0.0, 0.50000000000, 0.50000000000, 0.00003051758};

static MP2DEC_SIMD_ALIGN FixedReal hcos_64[16],
     hcos_32[ 8],
     hcos_16[ 4],
     hcos_8 [ 2],
     hcos_4;

static MP2DEC_SIMD_ALIGN FixedReal hcos_8_11nn[4];
static MP2DEC_SIMD_ALIGN FixedReal hcos_4_1ng1ng[4];
static MP2DEC_SIMD_ALIGN FixedReal hcos_32_4_reverse[4];

#define DUMP_DEBUG(ch)	\
	printf(#ch ": ");\
	printf("%f %f %f %f\n", ch##0 * 10000.0f, ch##1 * 10000.0f, ch##2 * 10000.0f, ch##3 * 10000.0f);\
	printf("   %f %f %f %f\n", ch##4 * 10000.0f, ch##5 * 10000.0f, ch##6 * 10000.0f, ch##7 * 10000.0f);\
	printf("   %f %f %f %f\n", ch##8 * 10000.0f, ch##9 * 10000.0f, ch##a * 10000.0f, ch##b * 10000.0f);\
	printf("   %f %f %f %f\n", ch##c * 10000.0f, ch##d * 10000.0f, ch##e * 10000.0f, ch##f * 10000.0f)

static void computebuffer(FixedReal *fraction, FixedReal buffer[2][MAXHISTORY][MAXSUBBAND], int currentcalcbuffer, int calcbufferoffset)
{
	FixedReal p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,pa,pb,pc,pd,pe,pf;
	FixedReal q0,q1,q2,q3,q4,q5,q6,q7,q8,q9,qa,qb,qc,qd,qe,qf;
	FixedReal *out1,*out2;

	out1 = buffer[currentcalcbuffer][calcbufferoffset];
	out2 = buffer[currentcalcbuffer^1][calcbufferoffset];
#define OUT1(v,t) out1[32-(v)]   =(-(out1[(v)]=t))
#define OUT2(v)   out2[32-(v)]=out2[(v)]

  // compute new values via a fast cosine transform:
  {
    register FixedReal *x=fraction;

    p0=x[ 0]+x[31];p1=x[ 1]+x[30];p2=x[ 2]+x[29];p3=x[ 3]+x[28];
    p4=x[ 4]+x[27];p5=x[ 5]+x[26];p6=x[ 6]+x[25];p7=x[ 7]+x[24];
    p8=x[ 8]+x[23];p9=x[ 9]+x[22];pa=x[10]+x[21];pb=x[11]+x[20];
    pc=x[12]+x[19];pd=x[13]+x[18];pe=x[14]+x[17];pf=x[15]+x[16];
  }

  q0=p0+pf;q1=p1+pe;q2=p2+pd;q3=p3+pc;
  q4=p4+pb;q5=p5+pa;q6=p6+p9;q7=p7+p8;
  q8=hcos_32[0]*(p0-pf);q9=hcos_32[1]*(p1-pe);
  qa=hcos_32[2]*(p2-pd);qb=hcos_32[3]*(p3-pc);
  qc=hcos_32[4]*(p4-pb);qd=hcos_32[5]*(p5-pa);
  qe=hcos_32[6]*(p6-p9);qf=hcos_32[7]*(p7-p8);

  p0=q0+q7;p1=q1+q6;p2=q2+q5;p3=q3+q4;
  p4=hcos_16[0]*(q0-q7);p5=hcos_16[1]*(q1-q6);
  p6=hcos_16[2]*(q2-q5);p7=hcos_16[3]*(q3-q4);
  p8=q8+qf;p9=q9+qe;pa=qa+qd;pb=qb+qc;
  pc=hcos_16[0]*(q8-qf);pd=hcos_16[1]*(q9-qe);
  pe=hcos_16[2]*(qa-qd);pf=hcos_16[3]*(qb-qc);
  
  q0=p0+p3;q1=p1+p2;q2=hcos_8[0]*(p0-p3);q3=hcos_8[1]*(p1-p2);
  q4=p4+p7;q5=p5+p6;q6=hcos_8[0]*(p4-p7);q7=hcos_8[1]*(p5-p6);
  q8=p8+pb;q9=p9+pa;qa=hcos_8[0]*(p8-pb);qb=hcos_8[1]*(p9-pa);
  qc=pc+pf;qd=pd+pe;qe=hcos_8[0]*(pc-pf);qf=hcos_8[1]*(pd-pe);

  p0=q0+q1;p1=hcos_4*(q0-q1);p2=q2+q3;p3=hcos_4*(q2-q3);
  p4=q4+q5;p5=hcos_4*(q4-q5);p6=q6+q7;p7=hcos_4*(q6-q7);
  p8=q8+q9;p9=hcos_4*(q8-q9);pa=qa+qb;pb=hcos_4*(qa-qb);
  pc=qc+qd;pd=hcos_4*(qc-qd);pe=qe+qf;pf=hcos_4*(qe-qf);
  
  {
    register FixedReal tmp;

    tmp=p6+p7;
    OUT2(4)=-(p5+tmp);
    OUT2(12)=-(p4+tmp);
    tmp=pb+pf;
    OUT1(10,tmp);
    OUT1(6,pd+tmp);
    tmp=pe+pf;
    OUT2(14)=-(p8+pc+tmp);
    OUT2(2)=-(p9+pd+tmp);
    tmp+=pa+pb;
    OUT2(6)=-(pd+tmp);
    OUT2(10)=-(pc+tmp);
    OUT1(2,p9+pd+pf);
    OUT1(4,p5+p7);
    OUT2(16)=-p0;
    out2[0]=-(out1[0]=p1);
    OUT1( 8,p3);
    OUT1(12,p7);
    OUT1(14,pf);
    OUT2(8)=-(p2+p3);
  }

  {
    register FixedReal *x=fraction;

    p0=hcos_64[ 0]*(x[ 0]-x[31]);p1=hcos_64[ 1]*(x[ 1]-x[30]);
    p2=hcos_64[ 2]*(x[ 2]-x[29]);p3=hcos_64[ 3]*(x[ 3]-x[28]);
    p4=hcos_64[ 4]*(x[ 4]-x[27]);p5=hcos_64[ 5]*(x[ 5]-x[26]);
    p6=hcos_64[ 6]*(x[ 6]-x[25]);p7=hcos_64[ 7]*(x[ 7]-x[24]);
    p8=hcos_64[ 8]*(x[ 8]-x[23]);p9=hcos_64[ 9]*(x[ 9]-x[22]);
    pa=hcos_64[10]*(x[10]-x[21]);pb=hcos_64[11]*(x[11]-x[20]);
    pc=hcos_64[12]*(x[12]-x[19]);pd=hcos_64[13]*(x[13]-x[18]);
    pe=hcos_64[14]*(x[14]-x[17]);pf=hcos_64[15]*(x[15]-x[16]);
  }
  
  q0=p0+pf;q1=p1+pe;q2=p2+pd;q3=p3+pc;
  q4=p4+pb;q5=p5+pa;q6=p6+p9;q7=p7+p8;
  q8=hcos_32[0]*(p0-pf);q9=hcos_32[1]*(p1-pe);
  qa=hcos_32[2]*(p2-pd);qb=hcos_32[3]*(p3-pc);
  qc=hcos_32[4]*(p4-pb);qd=hcos_32[5]*(p5-pa);
  qe=hcos_32[6]*(p6-p9);qf=hcos_32[7]*(p7-p8);
  
  p0=q0+q7;p1=q1+q6;p2=q2+q5;p3=q3+q4;
  p4=hcos_16[0]*(q0-q7);p5=hcos_16[1]*(q1-q6);
  p6=hcos_16[2]*(q2-q5);p7=hcos_16[3]*(q3-q4);
  p8=q8+qf;p9=q9+qe;pa=qa+qd;pb=qb+qc;
  pc=hcos_16[0]*(q8-qf);pd=hcos_16[1]*(q9-qe);
  pe=hcos_16[2]*(qa-qd);pf=hcos_16[3]*(qb-qc);
  
  q0=p0+p3;q1=p1+p2;q2=hcos_8[0]*(p0-p3);q3=hcos_8[1]*(p1-p2);
  q4=p4+p7;q5=p5+p6;q6=hcos_8[0]*(p4-p7);q7=hcos_8[1]*(p5-p6);
  q8=p8+pb;q9=p9+pa;qa=hcos_8[0]*(p8-pb);qb=hcos_8[1]*(p9-pa);
  qc=pc+pf;qd=pd+pe;qe=hcos_8[0]*(pc-pf);qf=hcos_8[1]*(pd-pe);
  
  p0=q0+q1;p1=hcos_4*(q0-q1);
  p2=q2+q3;p3=hcos_4*(q2-q3);
  p4=q4+q5;p5=hcos_4*(q4-q5);
  p6=q6+q7;p7=hcos_4*(q6-q7);
  p8=q8+q9;p9=hcos_4*(q8-q9);
  pa=qa+qb;pb=hcos_4*(qa-qb);
  pc=qc+qd;pd=hcos_4*(qc-qd);
  pe=qe+qf;pf=hcos_4*(qe-qf);
  
  {
    FixedReal tmp;

    tmp=pd+pf;
    OUT1(5,p5+p7+pb+tmp);
    tmp+=p9;
    OUT1(1,p1+tmp);
    OUT2(1)=-(p1+pe+tmp);
    tmp+=p5+p7;
    OUT1(3,tmp);
    OUT2(3)=-(p6+pe+tmp);
    tmp=pa+pb+pc+pd+pe+pf;
    OUT2(7)=-(p2+p3+tmp-pc);
    OUT2(11)=-(p4+p6+p7+tmp-pd);
    OUT2(5)=-(p5+p6+p7+tmp-pc);
    OUT2(9)=-(p2+p3+tmp-pd);
    tmp=p8+pc+pe+pf;
    OUT2(15)=-(p0+tmp);
    OUT2(13)=-(p4+p6+p7+tmp);
    tmp=pb+pf;
    OUT1(11,p7+tmp);
    tmp+=p3;
    OUT1( 9,tmp);
    OUT1( 7,pd+tmp);
    OUT1(13,p7+pf);
    OUT1(15,pf);
  }
}

inline __m128 sse_reverse_ps(__m128 v)
{
	return _mm_shuffle_ps(v, v, 0x1b);
}

inline __m128 sse_combine_lo(__m128 v)
{
	return _mm_movelh_ps(v, v);
}

inline __m128 sse_combine_and_reverse_hi(__m128 v)
{
	return _mm_xor_ps(_mm_set_ps(-0.0f, -0.0f, 0.0f, 0.0f), _mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 3, 2, 3)));
}

inline __m128 sse_reorder_1n03n2(__m128 v)
{
	return _mm_xor_ps(_mm_set_ps(-0.0f, 0.0f, -0.0f, 0.0f), _mm_shuffle_ps(v, v, _MM_SHUFFLE(2, 3, 0, 1)));
}
#define SSE_UNPACK(ch, n0,n1,n2,n3)	\
	FixedReal ch##n0 = ch##n0##n1##n2##n3.m128_f32[0];\
	FixedReal ch##n1 = ch##n0##n1##n2##n3.m128_f32[1];\
	FixedReal ch##n2 = ch##n0##n1##n2##n3.m128_f32[2];\
	FixedReal ch##n3 = ch##n0##n1##n2##n3.m128_f32[3]

#define SSE_UNPACK_ALL(ch)	\
	SSE_UNPACK(ch,0,1,2,3);\
	SSE_UNPACK(ch,4,5,6,7);\
	SSE_UNPACK(ch,8,9,a,b);\
	SSE_UNPACK(ch,c,d,e,f)

#define computebuffer_ct \
do {\
  __m128 hcos_32_0 = _mm_load_ps(hcos_32 + 0);\
  __m128 l_hcos_32_4_reverse = _mm_load_ps(hcos_32_4_reverse);\
  __m128 q0123 = _mm_add_ps(p0123, sse_reverse_ps(pcdef));\
  __m128 q_rev_4567 = _mm_add_ps(sse_reverse_ps(p4567), p89ab);\
  __m128 q89ab = _mm_mul_ps(hcos_32_0, _mm_sub_ps(p0123, sse_reverse_ps(pcdef)));\
  __m128 q_rev_cdef = _mm_mul_ps(l_hcos_32_4_reverse, _mm_sub_ps(sse_reverse_ps(p4567), p89ab));\
  \
  __m128 hcos_16_0 = _mm_load_ps(hcos_16);\
  p0123 = _mm_add_ps(q0123, q_rev_4567);\
  p4567 = _mm_mul_ps(hcos_16_0, _mm_sub_ps(q0123, q_rev_4567));\
  p89ab = _mm_add_ps(q89ab, q_rev_cdef);\
  pcdef = _mm_mul_ps(hcos_16_0, _mm_sub_ps(q89ab, q_rev_cdef));\
  \
  __m128 hcos8_11nn_0 = _mm_load_ps(hcos_8_11nn);\
  q0123 = _mm_mul_ps(hcos8_11nn_0, _mm_add_ps(sse_combine_lo(p0123), sse_combine_and_reverse_hi(p0123)));\
  __m128 q4567 = _mm_mul_ps(hcos8_11nn_0, _mm_add_ps(sse_combine_lo(p4567), sse_combine_and_reverse_hi(p4567)));\
  q89ab = _mm_mul_ps(hcos8_11nn_0, _mm_add_ps(sse_combine_lo(p89ab), sse_combine_and_reverse_hi(p89ab)));\
  __m128 qcdef = _mm_mul_ps(hcos8_11nn_0, _mm_add_ps(sse_combine_lo(pcdef), sse_combine_and_reverse_hi(pcdef)));\
  \
  __m128 l_hcos4_1ng1ng = _mm_load_ps(hcos_4_1ng1ng);\
  p0123 = _mm_mul_ps(l_hcos4_1ng1ng, _mm_add_ps(q0123, sse_reorder_1n03n2(q0123)));\
  p4567 = _mm_mul_ps(l_hcos4_1ng1ng, _mm_add_ps(q4567, sse_reorder_1n03n2(q4567)));\
  p89ab = _mm_mul_ps(l_hcos4_1ng1ng, _mm_add_ps(q89ab, sse_reorder_1n03n2(q89ab)));\
  pcdef = _mm_mul_ps(l_hcos4_1ng1ng, _mm_add_ps(qcdef, sse_reorder_1n03n2(qcdef)));\
} while(0)

static void computebuffer_sse(FixedReal *fraction, FixedReal buffer[2][MAXHISTORY][MAXSUBBAND], int currentcalcbuffer, int calcbufferoffset)
{
	
	__m128 frac0123 = _mm_load_ps(fraction+ 0);
	__m128 frac4567 = _mm_load_ps(fraction+ 4);
	__m128 frac89ab = _mm_load_ps(fraction+ 8);
	__m128 fraccdef = _mm_load_ps(fraction+12);
	__m128 frac0123_hi_rev = sse_reverse_ps(_mm_load_ps(fraction+16));
	__m128 frac4567_hi_rev = sse_reverse_ps(_mm_load_ps(fraction+20));
	__m128 frac89ab_hi_rev = sse_reverse_ps(_mm_load_ps(fraction+24));
	__m128 fraccdef_hi_rev = sse_reverse_ps(_mm_load_ps(fraction+28));

	__m128 p0123;
	__m128 p4567;
	__m128 p89ab;
	__m128 pcdef;

	FixedReal *out1,*out2;

	out1 = buffer[currentcalcbuffer][calcbufferoffset];
	out2 = buffer[currentcalcbuffer^1][calcbufferoffset];

  // compute new values via a fast cosine transform:
  {
    const FixedReal *x=fraction;

	p0123 = _mm_add_ps(frac0123, fraccdef_hi_rev);
	p4567 = _mm_add_ps(frac4567, frac89ab_hi_rev);
	p89ab = _mm_add_ps(frac89ab, frac4567_hi_rev);
	pcdef = _mm_add_ps(fraccdef, frac0123_hi_rev);
  }

  computebuffer_ct;
  
  {
    SSE_UNPACK_ALL(p);
	
    FixedReal tmp;

    tmp=p6+p7;
    OUT2(4)=-(p5+tmp);
    OUT2(12)=-(p4+tmp);
    tmp=pb+pf;
    OUT1(10,tmp);
    OUT1(6,pd+tmp);
    tmp=pe+pf;
    OUT2(14)=-(p8+pc+tmp);
    OUT2(2)=-(p9+pd+tmp);
    tmp+=pa+pb;
    OUT2(6)=-(pd+tmp);
    OUT2(10)=-(pc+tmp);
    OUT1(2,p9+pd+pf);
    OUT1(4,p5+p7);
    OUT2(16)=-p0;
    out2[0]=-(out1[0]=p1);
    OUT1( 8,p3);
    OUT1(12,p7);
    OUT1(14,pf);
    OUT2(8)=-(p2+p3);
  }

  {
    register FixedReal *x=fraction;

	p0123 = _mm_mul_ps(_mm_load_ps(hcos_64+ 0), _mm_sub_ps(frac0123, fraccdef_hi_rev));
	p4567 = _mm_mul_ps(_mm_load_ps(hcos_64+ 4), _mm_sub_ps(frac4567, frac89ab_hi_rev));
	p89ab = _mm_mul_ps(_mm_load_ps(hcos_64+ 8), _mm_sub_ps(frac89ab, frac4567_hi_rev));
	pcdef = _mm_mul_ps(_mm_load_ps(hcos_64+12), _mm_sub_ps(fraccdef, frac0123_hi_rev));
  }

  computebuffer_ct;

  {
    SSE_UNPACK_ALL(p);

    float tmp=pd+pf;

    OUT1(5,p5+p7+pb+tmp);
    tmp+=p9;
    OUT1(1,p1+tmp);
    OUT2(1)=-(p1+pe+tmp);
    tmp+=p5+p7;
    OUT1(3,tmp);
    OUT2(3)=-(p6+pe+tmp);
    tmp=pa+pb+pc+pd+pe+pf;
    OUT2(7)=-(p2+p3+tmp-pc);
    OUT2(11)=-(p4+p6+p7+tmp-pd);
    OUT2(5)=-(p5+p6+p7+tmp-pc);
    OUT2(9)=-(p2+p3+tmp-pd);
    tmp=p8+pc+pe+pf;
    OUT2(15)=-(p0+tmp);
    OUT2(13)=-(p4+p6+p7+tmp);
    tmp=pb+pf;
    OUT1(11,p7+tmp);
    tmp+=p3;
    OUT1( 9,tmp);
    OUT1( 7,pd+tmp);
    OUT1(13,p7+pf);
    OUT1(15,pf);
  }
}


#define SAVE \
		raw=(int)(MP2DEC_MAR(r, scalefactor)); \
		if(raw>MAXSCALE)raw=MAXSCALE;else if(raw<MINSCALE)raw=MINSCALE; \
		**output = MP2DEC_STATIC_CAST(mp2dec_sint16, raw);\
		(*output)++;
#define OS  r=*vp * *dp++
#define XX  vp+=15;r+=*vp * *dp++
#define OP  r+=*--vp * *dp++

static void generate(int currentcalcbuffer, int calcbufferoffset, FixedReal calcbuffer[2][MAXHISTORY][MAXSUBBAND], int scalefactor,
	mp2dec_sint16 **outputp)
{
  FixedReal samples[MAXSUBBAND];

  for(int band=0;band<MAXSUBBAND;band++)
	  samples[band] = FixedReal(0.0);

  const FixedReal (*currentcb)[32] = calcbuffer[currentcalcbuffer];

  for(int hist=0;hist<16;hist++)
  {
    const FixedReal *filterp = filter[hist];
    const FixedReal *cbp = currentcb[(calcbufferoffset - hist) & 15];

    for(int sample=0;sample<MAXSUBBAND;sample++)
      samples[sample] += filterp[sample] * cbp[sample];
  }

  mp2dec_sint16 *output = *outputp;
  for(int sample=0;sample<MAXSUBBAND;sample++)
  {
    int raw=(int)(MP2DEC_MAR(samples[sample], scalefactor));
    if(raw>MAXSCALE)
      raw=MAXSCALE;
    else if(raw<MINSCALE)
      raw=MINSCALE;
    *output++ = MP2DEC_STATIC_CAST(mp2dec_sint16, raw);
  }
  *outputp = output;
}

#ifdef MP2DEC_ENABLE_SSE2


#define SSE2_FAST_TRANSPOSE(r0, r1, r2, r3)	\
	{\
		__m128 tmp0 = _mm_unpacklo_ps(r0, r1);\
		__m128 tmp1 = _mm_unpacklo_ps(r2, r3);\
		__m128 tmp2 = _mm_unpackhi_ps(r0, r1);\
		__m128 tmp3 = _mm_unpackhi_ps(r2, r3);\
		r0 = _mm_movelh_ps(tmp0, tmp1);\
		r1 = _mm_movehl_ps(tmp1, tmp0);\
		r2 = _mm_movelh_ps(tmp2, tmp3);\
		r3 = _mm_movehl_ps(tmp3, tmp2);\
	}


static void generate_sse2(int currentcalcbuffer, int calcbufferoffset, FixedReal calcbuffer[2][MAXHISTORY][MAXSUBBAND], int scalefactor,
	mp2dec_sint16 **outputp)
{
  __m128 sampleCluster[MAXSUBBAND/4];

  const FixedReal (*currentcb)[32] = calcbuffer[currentcalcbuffer];

  {
    const FixedReal *filterp = filter[0];
    const FixedReal *cbp = currentcb[calcbufferoffset];

    for(int sample=0;sample<MAXSUBBAND/4;sample++)
      sampleCluster[sample] = _mm_mul_ps(_mm_load_ps(filterp + sample*4), _mm_load_ps(cbp + sample*4));
  }

  for(int hist=1;hist<16;hist++)
  {
    const FixedReal *filterp = filter[hist];
    const FixedReal *cbp = currentcb[(calcbufferoffset - hist) & 15];

    for(int sample=0;sample<MAXSUBBAND/4;sample++)
      sampleCluster[sample] = _mm_add_ps(sampleCluster[sample], _mm_mul_ps(_mm_load_ps(filterp + sample*4), _mm_load_ps(cbp + sample*4)));
  }

  mp2dec_sint16 *output = *outputp;

  // Rescale
  for(int sample=0;sample<MAXSUBBAND/8;sample++)
  {
    __m128i ipack0 = _mm_cvtps_epi32(sampleCluster[sample*2+0]);
    __m128i ipack1 = _mm_cvtps_epi32(sampleCluster[sample*2+1]);
    __m128i packed = _mm_packs_epi32(ipack0, ipack1);
    _mm_store_si128(reinterpret_cast<__m128i*>(output), packed);
	output += 8;
  }
  *outputp = output;
}
#endif


static void subbandsynthesis(FixedReal *fraction, int *currentcalcbuffer, int *calcbufferoffset, FixedReal calcbuffer[2][MAXHISTORY][MAXSUBBAND], mp2dec_sint16 **pOutput)
{
#ifdef MP2DEC_ENABLE_SSE2
	computebuffer_sse(fraction, calcbuffer, *currentcalcbuffer, *calcbufferoffset);
	generate_sse2(*currentcalcbuffer, *calcbufferoffset, calcbuffer, SCALE, pOutput);
#else
	computebuffer(fraction, calcbuffer, *currentcalcbuffer, *calcbufferoffset, false);
	generate(*currentcalcbuffer, *calcbufferoffset, calcbuffer, SCALE, pOutput);
#endif

	if(*calcbufferoffset < 15)
		(*calcbufferoffset)++;
	else
		*calcbufferoffset = 0;

	(*currentcalcbuffer) ^=1;
}

// For header
enum
{
	mp2dec_version_mpeg2,
	mp2dec_version_mpeg1
};

enum
{
	mp2dec_mode_fullstereo,
	mp2dec_mode_joint,
	mp2dec_mode_dual,
	mp2dec_mode_single
};

static const int tbl_frequencies[2][3]=
{
	{22050,24000,16000}, // MPEG 2
	{44100,48000,32000}  // MPEG 1
};

static const int tbl_bitrates[2][3][15]=
{
	// MPEG 2
	{{0,32,48,56,64,80,96,112,128,144,160,176,192,224,256},
	 {0,8,16,24,32,40,48,56,64,80,96,112,128,144,160},
	 {0,8,16,24,32,40,48,56,64,80,96,112,128,144,160}},

	// MPEG 1
	{{0,32,64,96,128,160,192,224,256,288,320,352,384,416,448},
	 {0,32,48,56,64,80,96,112,128,160,192,224,256,320,384},
	 {0,32,40,48,56,64,80,96,112,128,160,192,224,256,320}}
};

static const int tbl_bitalloclength[2][MAXSUBBAND]=
{{4,4,3,3,3,3,3,3,3,3,3,3,3,3,3,3, 3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3},
 {4,4,4,4,4,4,4,4,4,4,4,3,3,3,3,3, 3,3,3,3,3,3,3,2,2,2,2,2,2,2,2,2}};

static struct mp2dec_auto_init
{
	mp2dec_auto_init(void)
	{
		const FixedReal group5base[3] = { -2.0/3.0, 0.0, 2.0/3.0 };
		const FixedReal group7base[5] = { -0.8, -0.4, 0.0, 0.4, 0.8 };
		const FixedReal group10base[9] = { -8.0/9.0, -6.0/9.0, -4.0/9.0, -2.0/9.0, 0.0, 2.0/9.0, 4.0/9.0, 6.0/9.0, 8.0/9.0 };

		int i;

		for(i=0;i<16;i++) hcos_64[i] = (FixedReal)
				(1.0/(2.0*cos(M_PI*MP2DEC_STATIC_CAST(double, i*2+1)/64.0)));
		for(i=0;i< 8;i++) hcos_32[i] = (FixedReal)
				(1.0/(2.0*cos(M_PI*MP2DEC_STATIC_CAST(double, i*2+1)/32.0)));
		for(i=0;i< 4;i++) hcos_16[i] = (FixedReal)
				(1.0/(2.0*cos(M_PI*MP2DEC_STATIC_CAST(double, i*2+1)/16.0)));
		for(i=0;i< 2;i++) hcos_8 [i] = (FixedReal)
				(1.0/(2.0*cos(M_PI*MP2DEC_STATIC_CAST(double, i*2+1)/ 8.0)));
		hcos_4 = (FixedReal)(1.0f / (2.0f * cos( M_PI * 1.0 / 4.0 )));

		hcos_8_11nn[0] = 1.0f;
		hcos_8_11nn[1] = 1.0f;
		hcos_8_11nn[2] = hcos_8[0];
		hcos_8_11nn[3] = hcos_8[1];

		hcos_4_1ng1ng[0] = 1.0f;
		hcos_4_1ng1ng[1] = -hcos_4;
		hcos_4_1ng1ng[2] = 1.0f;
		hcos_4_1ng1ng[3] = -hcos_4;

		for(i=0;i< 4;i++)
			hcos_32_4_reverse[i] = hcos_32[7-i];

#ifdef MP2DEC_ENABLE_SSE2
		for(int hist=0;hist<MAXHISTORY;hist++)
			for(int band=0;band<MAXSUBBAND;band++)
				filter[hist][band] *= static_cast<FixedReal>(SCALE);
#endif

		for(i=0;i<64;i++)
			scalefactorstable[i] = FixedReal(2.0 * pow(0.5, double(i) / 3.0));


		for(i=0;i<27;i++)
		{
			group5bits[i*3+0] = group5base[i%3];
			group5bits[i*3+1] = group5base[(i/3)%3];
			group5bits[i*3+2] = group5base[i/9];
		}
		for(i=0;i<125;i++)
		{
			group7bits[i*3+0] = group7base[i%5];
			group7bits[i*3+1] = group7base[(i/5)%5];
			group7bits[i*3+2] = group7base[i/25];
		}
		for(i=0;i<729;i++)
		{
			group10bits[i*3+0] = group10base[i%9];
			group10bits[i*3+1] = group10base[(i/9)%9];
			group10bits[i*3+2] = group10base[i/81];
		}
	}
} mp2dec_auto_init_instance;

void mp2dec_bitstream_init(mp2dec_bitstream *bitstream, const void *data)
{
	bitstream->bytes = MP2DEC_STATIC_CAST(const mp2dec_uint8 *, data);
	bitstream->bitOffs = 0;
	bitstream->totalDec = 0;
}

mp2dec_uint mp2dec_bitstream_read(mp2dec_bitstream *bitstream, mp2dec_uint bits)
{
	mp2dec_uint result = 0;
	mp2dec_uint resultBitOffs = bits;
	mp2dec_uint streamBitOffs = bitstream->bitOffs;
	mp2dec_uint availBits;
	const mp2dec_uint8 *bytes = bitstream->bytes;
	mp2dec_uint originalBits = bits;	// Compiler optimization assist

	bitstream->totalDec += bits;

	if(bits == 1)
	{
		result = MP2DEC_STREAM_BIT(bytes[0], streamBitOffs);
		if(streamBitOffs == 7)
		{
			bitstream->bytes = bytes + 1;
			bitstream->bitOffs = 0;
		}
		else
			bitstream->bitOffs++;
		return result;
	}

	if(streamBitOffs != 0)
	{
		availBits = 8 - streamBitOffs;
		if(bits >= availBits)
		{
			// Longer than the first byte
			resultBitOffs -= availBits;
			result = (bytes[0] & MP2DEC_BIT_MASK(0, availBits)) << resultBitOffs;
			bits -= availBits;
			streamBitOffs = 0;
			bytes++;
		}
		else
		{
			result = (bytes[0] & MP2DEC_BIT_MASK((availBits - bits), bits)) >> (availBits - bits);
			streamBitOffs += bits;

			bitstream->bitOffs = streamBitOffs;
			return result;
		}
	}

	if(originalBits >= 8)	// Compiler optimization assist
	{
		while(bits >= 8)
		{
			resultBitOffs -= 8;
			result |= (bytes[0] << resultBitOffs);
			bits -= 8;
			bytes++;
		}
	}

	if(bits != 0)
	{
		result |= (bytes[0] & MP2DEC_BIT_MASK(8 - bits, bits)) >> (8 - bits);
		streamBitOffs = bits;
	}

	bitstream->bytes = bytes;
	bitstream->bitOffs = streamBitOffs;

	return result;
}

void mp2dec_init_decode_state(mp2dec_decode_properties *decodeProps)
{
	int i;
	memset(decodeProps->calcBuffersUnaligned, 0, sizeof(decodeProps->calcBuffersUnaligned));
	for(i=0;i<MAXCHANNEL;i++)
	{
		decodeProps->currentcalcbuffers[i] = 0;
		decodeProps->calcbufferoffsets[i] = 15;
	}
	decodeProps->cbAlignmentOffset = 0;
}

mp2dec_bool mp2dec_parse_header(const void *inBytes, mp2dec_header *outHeader, mp2dec_decode_properties *decodeProps)
{
	mp2dec_bitstream bitstream;

	mp2dec_bool inputStereo;
	mp2dec_uint channelBitrate;
	mp2dec_uint totalBitrate;
	mp2dec_uint tableIndex;
	mp2dec_uint numSubBands;
	mp2dec_uint sampleRate;
	mp2dec_uint stereoBound;

	mp2dec_bitstream_init(&bitstream, inBytes);

	outHeader->syncword = mp2dec_bitstream_read(&bitstream, 12);
	outHeader->version = mp2dec_bitstream_read(&bitstream, 1);
	outHeader->layer = mp2dec_bitstream_read(&bitstream, 2);
	outHeader->disableCRCFlag = mp2dec_bitstream_read(&bitstream, 1);
	outHeader->bitrateIndex = mp2dec_bitstream_read(&bitstream, 4);
	outHeader->samplerateIndex = mp2dec_bitstream_read(&bitstream, 2);
	outHeader->padding = mp2dec_bitstream_read(&bitstream, 1);
	outHeader->privateFlag = mp2dec_bitstream_read(&bitstream, 1);
	outHeader->mode = mp2dec_bitstream_read(&bitstream, 2);
	outHeader->modeExt = mp2dec_bitstream_read(&bitstream, 2);
	outHeader->copyrightFlag = mp2dec_bitstream_read(&bitstream, 1);
	outHeader->originalFlag = mp2dec_bitstream_read(&bitstream, 1);
	outHeader->emphasis = mp2dec_bitstream_read(&bitstream, 2);


	inputStereo = (outHeader->mode == mp2dec_mode_single) ? MP2DEC_FALSE : MP2DEC_TRUE;

	if(outHeader->syncword != 0xfff
		|| outHeader->bitrateIndex == 15
		|| outHeader->samplerateIndex == 3
		|| outHeader->layer != 2)
		return MP2DEC_FALSE;

	totalBitrate = channelBitrate = tbl_bitrates[outHeader->version][outHeader->layer - 1][outHeader->bitrateIndex];
	if(inputStereo)
		channelBitrate /= 2;

	if(outHeader->bitrateIndex < 3)
		tableIndex = 0;
	else
		tableIndex = 1;

	sampleRate = tbl_frequencies[outHeader->version][outHeader->samplerateIndex];

	if(outHeader->version == mp2dec_version_mpeg1)
	{
		if(channelBitrate >= 56 &&
			(sampleRate == 48000 || channelBitrate <= 80))
			numSubBands = 27;
		else if(sampleRate != 48000 && channelBitrate >= 96)
			numSubBands = 30;
		else if (sampleRate != 32000 && channelBitrate <= 48)
			numSubBands = 8;
		else
			numSubBands = 12;
	}
	else
		numSubBands = 30;

	if(outHeader->mode == mp2dec_mode_single)
		stereoBound = 0;
	else if(outHeader->mode == mp2dec_mode_joint)
		stereoBound = (outHeader->modeExt + 1) << 2;
	else
		stereoBound = numSubBands;

	if(stereoBound > numSubBands)
		stereoBound = numSubBands;

	// 1152 * 1000 / 8
	// totalBitrate * 1000 * 1152/ sampleRate
	decodeProps->frameBytes = (totalBitrate * 144000 / sampleRate) - 4;
	decodeProps->tableIndex = tableIndex;
	decodeProps->numSubBands = numSubBands;
	decodeProps->stereoBound = stereoBound;
	decodeProps->inputStereo = inputStereo;
	decodeProps->sampleRate = sampleRate;

	if(outHeader->padding)
		decodeProps->frameBytes++;

	return MP2DEC_TRUE;
}

template<int numChannels>
inline void mp2dec_decode_frame_ch(const mp2dec_header *header, mp2dec_decode_properties *decodeProps,
	FixedReal (*calcbuffers)[2][16][32], const void *frameData, mp2dec_sint16 output[2][1152])
{
	char fractionBuffer[sizeof(FixedReal)*MAXCHANNEL*3*MAXSUBBAND+MP2DEC_SIMD_SIZE];

	char *fractionAlignable = fractionBuffer + MP2DEC_ALIGN_OFFSET(fractionBuffer);
	
	FixedReal (*fraction)[3][MAXSUBBAND] = reinterpret_cast<FixedReal (*)[3][MAXSUBBAND]>(fractionAlignable);
	

	unsigned int bitalloc[MAXCHANNEL][MAXSUBBAND],
		scaleselector[MAXCHANNEL][MAXSUBBAND];
	FixedReal scalefactor[2][3][MAXSUBBAND];

	const FixedReal *group[MAXCHANNEL][MAXSUBBAND];
	unsigned int codelength[MAXCHANNEL][MAXSUBBAND];
	FixedReal factor[MAXCHANNEL][MAXSUBBAND];
	FixedReal c[MAXCHANNEL][MAXSUBBAND],d[MAXCHANNEL][MAXSUBBAND];

	mp2dec_sint16 *outputPointers[MAXCHANNEL];

	mp2dec_uint tableIndex = decodeProps->tableIndex;
	mp2dec_uint numSubBands = decodeProps->numSubBands;
	mp2dec_uint stereoBound = decodeProps->stereoBound;

	mp2dec_bitstream bitstream;

	
	mp2dec_bitstream_init(&bitstream, frameData);

	{
		int i;
		for(i=0;i<MAXCHANNEL;i++)
			outputPointers[i] = output[i];
	}

	if(!header->disableCRCFlag)
		mp2dec_bitstream_read(&bitstream, 16);

	if(numChannels == 1)
		stereoBound = 0;	// Compiler assist

	// Bitalloc
	{
		register mp2dec_uint i;
		register int ch;
		register const int *t = tbl_bitalloclength[tableIndex];

		if(PVERBOSE) printf("\nAllocs:\n");
		for(i=0;i<stereoBound;i++,t++)
			for(ch=0;ch<numChannels;ch++)
			{
				bitalloc[ch][i] = mp2dec_bitstream_read(&bitstream, *t);
				if(PVERBOSE) printf("%i:%i ", bitalloc[ch][i], *t);
			}
		for(;i<numSubBands;i++,t++)
		{
			mp2dec_uint l = mp2dec_bitstream_read(&bitstream, *t);
			for(ch=0;ch<numChannels;ch++)
				bitalloc[ch][i] = l;
		}
	}


	// Scale selector
	{
		int ch;
		mp2dec_uint i;
		if(PVERBOSE) printf("\nDScale selector\n");
		for(i=0;i<numSubBands;i++)
			for(ch=0;ch<numChannels;ch++)
				if(bitalloc[ch][i])
				{
					scaleselector[ch][i] = mp2dec_bitstream_read(&bitstream, 2);
					if(PVERBOSE) printf("%i ", scaleselector[ch][i]);
				}
	}
	
	if(PVERBOSE) printf("\nDScale factors\n");
	// Scale index
	{
		for(mp2dec_uint i=0;i<numSubBands;i++)
		{
			int ch;
			for(ch=0;ch<numChannels;ch++)
			{
				register int ba = bitalloc[ch][i];

				if(ba)
				{
					if(!tableIndex)
					{
						group[ch][i]=grouptableA[ba];
						codelength[ch][i]=codelengthtableA[ba];
						factor[ch][i]=factortableA[ba];
						c[ch][i]=ctableA[ba];
						d[ch][i]=dtableA[ba];
					}
					else
					{
						if(i <= 2)
						{
							group[ch][i]=grouptableB1[ba];
							codelength[ch][i]=codelengthtableB1[ba];
							factor[ch][i]=factortableB1[ba];
							c[ch][i]=ctableB1[ba];
							d[ch][i]=dtableB1[ba];
						}
						else
						{
							group[ch][i] = grouptableB234[ba];
							if(i <= 10)
							{
								codelength[ch][i] = codelengthtableB2[ba];
								factor[ch][i] = factortableB2[ba];
								c[ch][i] = ctableB2[ba];
								d[ch][i] = dtableB2[ba];
							}
							else if(i <= 22)
							{
								codelength[ch][i]=codelengthtableB3[ba];
								factor[ch][i]=factortableB3[ba];
								c[ch][i]=ctableB3[ba];
								d[ch][i]=dtableB3[ba];
							}
							else
							{
								codelength[ch][i]=codelengthtableB4[ba];
								factor[ch][i]=factortableB4[ba];
								c[ch][i]=ctableB4[ba];
								d[ch][i]=dtableB4[ba];
							}
						}
					}

					if(PVERBOSE) printf("%i:", scaleselector[ch][i]);
					switch(scaleselector[ch][i])
					{
					case 0:
						{
							mp2dec_uint32 bits = mp2dec_bitstream_read(&bitstream, 18);
							scalefactor[ch][0][i] = scalefactorstable[(bits & MP2DEC_BIT_MASK(12, 6)) >> 12];
							scalefactor[ch][1][i] = scalefactorstable[(bits & MP2DEC_BIT_MASK(6, 6)) >> 6];
							scalefactor[ch][2][i] = scalefactorstable[(bits & MP2DEC_BIT_MASK(0, 6))];
							if(PVERBOSE) printf("%i %i %i ", ((bits & MP2DEC_BIT_MASK(12, 6)) >> 12), ((bits & MP2DEC_BIT_MASK(6, 6)) >> 6), (bits & MP2DEC_BIT_MASK(0, 6)));
						}
						break;
					case 1:
						{
							mp2dec_uint32 bits = mp2dec_bitstream_read(&bitstream, 12);
							scalefactor[ch][0][i]=
								scalefactor[ch][1][i]=scalefactorstable[(bits & MP2DEC_BIT_MASK(6, 6)) >> 6];
							scalefactor[ch][2][i]=scalefactorstable[(bits & MP2DEC_BIT_MASK(0, 6))];
							if(PVERBOSE) printf("%i %i ", ((bits & MP2DEC_BIT_MASK(6, 6)) >> 6), (bits & MP2DEC_BIT_MASK(0, 6)));
						}
						break;
					case 2:
						{
							mp2dec_uint bits = mp2dec_bitstream_read(&bitstream, 6);
							scalefactor[ch][0][i]=
							   scalefactor[ch][1][i]=
							   scalefactor[ch][2][i]=scalefactorstable[bits];
							if(PVERBOSE) printf("%i ", bits);
						}
						break;
					case 3:
						{
							mp2dec_uint32 bits = mp2dec_bitstream_read(&bitstream, 12);
							scalefactor[ch][0][i]=scalefactorstable[(bits & MP2DEC_BIT_MASK(6, 6)) >> 6];
							scalefactor[ch][1][i]=
								scalefactor[ch][2][i]=scalefactorstable[(bits & MP2DEC_BIT_MASK(0, 6))];
							if(PVERBOSE) printf("%i %i ", ((bits & MP2DEC_BIT_MASK(6, 6)) >> 6), (bits & MP2DEC_BIT_MASK(0, 6)));
						}
						break;
					}
				}
			}	// ch
		}	//i
	}

	if(PVERBOSE) printf("\nDSample:\n");
	// Read Sample
	{
		int block;

		for(block=0;block<SCALEBLOCK;block++)
		{
			// Read Sample
			for(mp2dec_uint i=0;i<numSubBands;i++)
			{
				int ch;
				for(ch=0;ch<numChannels;ch++)
				{
					if(i >= stereoBound && ch > 0)
					{
						int j;
						for(j=0;j<3;j++)
							fraction[ch][j][i] = fraction[0][j][i];
						continue;
					}

					if(bitalloc[ch][i])
					{
						if(group[ch][i])
						{
							register const FixedReal *s;
							int code = mp2dec_bitstream_read(&bitstream, codelength[ch][i]);
							if(PVERBOSE) printf("%i:%i ", codelength[ch][i], code, codelength[ch][i]);

							code += code << 1;
							if (code > 2184)
								code = 2184;
							s = group[ch][i] + code;

							fraction[ch][0][i] = s[0];
							fraction[ch][1][i] = s[1];
							fraction[ch][2][i] = s[2];
						}
						else
						{
							mp2dec_uint l = codelength[ch][i];
							mp2dec_uint32 bits = mp2dec_bitstream_read(&bitstream, l*3);
							fraction[ch][0][i] =
								MP2DEC_STATIC_CAST(FixedReal17, (bits & MP2DEC_BIT_MASK(l*2, l)) >> (l*2)) * factor[ch][i] - FixedReal17(1.0);
							fraction[ch][1][i] =
								MP2DEC_STATIC_CAST(FixedReal17, (bits & MP2DEC_BIT_MASK(l, l)) >> l) * factor[ch][i] - FixedReal17(1.0);
							fraction[ch][2][i] =
								MP2DEC_STATIC_CAST(FixedReal17, (bits & MP2DEC_BIT_MASK(0, l))) * factor[ch][i] - FixedReal17(1.0);

							if(PVERBOSE) printf("%i:%i ", (bits & MP2DEC_BIT_MASK(l*2, l)) >> (l*2), l);
							if(PVERBOSE) printf("%i:%i ", (bits & MP2DEC_BIT_MASK(l, l)) >> l, l);
							if(PVERBOSE) printf("%i:%i ", (bits & MP2DEC_BIT_MASK(0, l)), l);
						}
					}
					else
						fraction[ch][0][i] = fraction[ch][1][i] = fraction[ch][2][i] = 0.0;
				}
			}

			//Fraction
			{
				int ch;
				for(ch=0;ch<numChannels;ch++)
				{
					for(mp2dec_uint i=0;i<numSubBands;i++)
					{
						if(bitalloc[ch][i])
						{
							register FixedReal t = scalefactor[ch][block/4][i];

							if(!group[ch][i])
							{
								fraction[ch][0][i] = (fraction[ch][0][i] + d[ch][i]) * c[ch][i];
								fraction[ch][1][i] = (fraction[ch][1][i] + d[ch][i]) * c[ch][i];
								fraction[ch][2][i] = (fraction[ch][2][i] + d[ch][i]) * c[ch][i];
							}

							fraction[ch][0][i] *= t;
							fraction[ch][1][i] *= t;
							fraction[ch][2][i] *= t;
						}
					}
				}
			}

			{
				int i;
				for(i=numSubBands;i<MAXSUBBAND;i++)
				{
					int ch;
					for(ch=0;ch<numChannels;ch++)
						fraction[ch][0][i] = fraction[ch][1][i] = fraction[ch][2][i]=0.0;
				}
			}

			{
				int i;
				for(i=0;i<3;i++)
				{
					int ch;
					for(ch=0;ch<numChannels;ch++)
						subbandsynthesis(fraction[ch][i], &decodeProps->currentcalcbuffers[ch], &decodeProps->calcbufferoffsets[ch],
							calcbuffers[ch], &outputPointers[ch]);
				}
			}
		}
	}
}

void mp2dec_decode_frame(const mp2dec_header *header, mp2dec_decode_properties *decodeProps,
	const void *frameData, mp2dec_sint16 output[2][1152])
{
	unsigned int expectedOffset = static_cast<unsigned int>(MP2DEC_ALIGN_OFFSET(decodeProps->calcBuffersUnaligned));
	if(expectedOffset != decodeProps->cbAlignmentOffset)
	{
		memmove(decodeProps->calcBuffersUnaligned + expectedOffset,
			decodeProps->calcBuffersUnaligned + decodeProps->cbAlignmentOffset,
			sizeof(decodeProps->calcBuffersUnaligned) - MP2DEC_SIMD_SIZE);
		decodeProps->cbAlignmentOffset = expectedOffset;
	}

	FixedReal (*calcbuffers)[2][16][32] = reinterpret_cast<FixedReal (*)[2][16][32]>(decodeProps->calcBuffersUnaligned + expectedOffset);

	if(decodeProps->inputStereo)
		mp2dec_decode_frame_ch<2>(header, decodeProps, calcbuffers, frameData, output);
	else
		mp2dec_decode_frame_ch<1>(header, decodeProps, calcbuffers, frameData, output);
}

void mp2dec_mix(const mp2dec_sint16 input[2][1152], mp2dec_sint16 mixed[1152*2])
{
	int i,j;
	for(i=0;i<1152;i++)
		for(j=0;j<2;j++)
			mixed[i*2+j] = input[j][i];
}
